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We analyze the problem of eliminating finite-size errors from quantum Monte Carlo (QMC) energy 
data. We demonstrate that both (i) adding a recently proposed 1 finite-size correction to the Ewald 
■ energy and (ii) using the model periodic Coulomb (MPC) interaction^* 4 - are good solutions to the 

| problem of removing finite-size effects from the interaction energy in cubic systems, provided the 

exchange-correlation (XC) hole has converged with respect to system size. However, we find that the 
MPC interaction distorts the XC hole in finite systems, implying that the Ewald interaction should 
be used to generate the configuration distribution. The finite-size correction of Ref. [f| is shown to 
be incomplete in systems of low symmetry. Beyond-leading-order corrections to the kinetic energy 
in . are found to be necessary at intermediate and high densities, and we investigate the effect of adding 

such corrections to QMC data for the homogeneous electron gas. We analyze finite-size errors 
in two-dimensional systems and show that the leading-order behavior differs from that which has 
O ■ hitherto been supposed. We compare the efficiency of different twist- averaging methods for reducing 

single-particle finite-size errors and we examine the performance of various finite-size extrapolation 
formulas. Finally, we investigate the system-size scaling of biases in diffusion QMC. 

PACS numbers: 02.70.Ss, 71.15.Nc, 71.10.Ca 

g ■ I. INTRODUCTION 

Continuum quantum Monte Carlo (QMC) techniques^ enable the total energies of many-electron systems to be 
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calculated to very high accuracy. QMC simulations of condensed matter are usually performed using finite simulation 
cells subject to periodic boundary conditions. The energy per particle is calculated at several different system sizes 
and the results are extrapolated to infinite system size. Unfortunately, this process introduces errors into the QMC 
results. Indeed, for simple systems such as the homogeneous electron gas (HEG), finite-size extrapolation is believed 
to be the largest single source of error in QMC data. 

In this paper we address a number of outstanding problems associated with finite-size extrapolation. We discuss 
the physics of finite-size effects in Sec. HU In Sec. [TTT] we discuss the use of twist-averaged boundary conditions^ to 
reduce errors caused by momentum quantization in finite simulation cells. In Sec. IIVI we give results illustrating 
that recently proposed methods for correcting the Ewald energy 1 -^ are essentially equivalent to the use of the model 
, periodic Coulomb (MPC) interaction 2,3,4 in QMC simulations. In Sec. [V] we discuss various complications posed by 
low-symmetry systems. In Sec. IVII we demonstrate that the finite-size correction to the kinetic energy (KE) proposed 
in Ref. [l| is incomplete and that higher-order terms cannot be neglected at typical metallic densities. We analyze 
finite-size errors in 2D-periodic systems in Sec. lVIIl In Sec. IVIlJl wc investigate the performance of different finite-size 
extrapolation formulas. In Sec. IIXI we examine the size-dependence of biases in QMC energies. Finally we draw our 
• i— i , conclusions in Sec. IXl 

Hartree atomic units (a.u.) are used throughout, in which the Dirac constant, the magnitude of the electronic 
charge, the electronic mass, and An times the permittivity of free space are unity: K = \e\ = m e = 47reo = 1. All 
our QMC calculations were carried out using the CASINO code4 We have made use of the variational and diffusion 
quantum Monte Carlo (VMC and DMC) methods^ Throughout, we specify the density of a HEG by quoting the 
radius r s of the sphere (circle in 2D) that contains one electron on average. 
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II. PHYSICS OF FINITE-SIZE EFFECTS 



A. Components of the total energy 

The total energy of a periodic many-electron system can be divided into (i) the KE, (ii) the electron-electron 
interaction energy, and (iii) the electron-ion interaction energy (we include the interaction of the electrons with any 
other external fields in this term). The electron-electron interaction energy may be subdivided into the Hartree and 
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exchange-correlation (XC) energies. Assuming the electrostatic potential to be periodic, the former is the Coulomb 
energy due to the periodic charge density, and the latter is the remainder of the electron-electron interaction energy, 
which arises from the correlation of electron motions and the antisymmetry of the many-electron wave function. The 
electron-ion interaction energy and the Hartree energy depend only on the electronic charge density, which has the 
periodicity of the primitive cell and is rapidly convergent with respect to system size; hence the finite-size errors in 
these energy components are generally negligible. By contrast, the finite-size errors in the XC energy and the KE can 
be very substantial. We analyze the physics underlying these errors in the rest of this section. 



B. Simulation and primitive unit cells for crystalline solids 



Suppose we wish to calculate the energy per particle of a periodic solid. In one-electron theories we can often reduce 
the problem to the primitive unit cell and integrate over the first Brillouin zone. Reduction to the primitive cell is 
not possible in many-body calculations because correlation effects may be long-ranged, and hence such calculations 
must be performed in simulation cells consisting of several primitive cells. Periodic boundary conditions are imposed 
across the simulation cell. 

Suppose the simulation cell contains N electrons and let {ri, . . . ,rjv} be the electron coordinates. The simulation- 
cell Hamiltonian H satisfies 

J?(ri,...,ri + R fl ,...,rjv) = H(r u . . . ,r u . . . ,r N ) Vi € {1, . . . , N} (1) 
iT(ri+R p ,...,r i + R p ,...,r iV + Rp) = £T(ri, . . . ,r i} . . . ,r N ), (2) 

where R s and R p are simulation-cell and primitive-cell lattice vectors. The first of these symmetries is an artifact 
of the periodicity imposed on the simulation cell. These translational symmetries lead to the many-body Bloch 
conditions 



* ks (ri,...,rjv) = Uk a (ri, . . . , r N ) exp ^ik s ■ ^ r^J (3) 
*kp(ri,...,rjv) = W / k p (ri,...,r^)exp ^ik p ■ I ' ( 4 ) 



where U has the periodicity of the simulation cell for every electron and W is invariant under the simultaneous 
translation of all electrons through R p i2il£ The use of a nonzero simulation-cell Bloch vector k s is sometimes described 
as the application of twisted boundary conditions^ The center-of-mass Bloch momentum k p may be restricted to the 
Brillouin zone corresponding to the primitive lattice, while the twist vector k s may be restricted to the smaller 
Brillouin zone corresponding to the simulation-cell lattice. From now on, we use G s and G p to denote vectors in the 
simulation-cell and primitive-cell reciprocal lattices, respectively. 



C. Single- particle finite-size errors 



In a finite simulation cell subject to periodic boundary conditions, each single-particle orbital can be taken to be of 
Bloch form V'k(r) = exp[ik • r]u k (r), where «k has the periodicity of the primitive cell and k lies on the grid of integer 
multiples of the simulation-cell reciprocal-lattice vectors within the first Brillouin zone of the primitive cell, the grid 
being offset from the origin by k s , so that k = k s + G s for some G s . Instead of integrating over single-particle orbitals 
inside the Fermi surface to calculate the HF KE and exchange energy, one therefore sums over a discrete set of k 
vectors when a finite cell is used. For metallic systems, the set of occupied ground-state orbitals depends on k s and 
hence calculated properties are nonanalytic functions of k s . As the system size is increased, the fineness of the grid 
of single-particle Bloch k vectors increases and the HF energy changes abruptly as shells of orbitals pass through the 
Fermi surface. 

Fluctuations in the QMC KE contain large "single-particle" contributions that are roughly proportional to the 
corresponding fluctuations in the HF KE. Hence HF energy data can be used to extrapolate QMC energies to infinite 
system size, as discussed in Sec. I Villi Note that a judicious choice of k s (e.g., the Baldereschi point— for insulators) 
can greatly reduce single-particle finite-size errors**™ The common choice of k s = generally maximizes shell-filling 
effects and is usually the worst possible value for estimating the total energy, although it does ensure that the wave 
function of the finite simulation cell can be chosen to have the full symmetry of the Hamiltonian. 



3 



D. Twist averaging 

Twist averaging within the canonical ensemble (CE) means taking the average of expectation values over all 
simulation-cell Bloch vectors k s in the first Brillouin zone of the simulation cell, i.e., over all offsets to the grid 
of k vectors, keeping the number of electrons fixedi^ 

At the HF level, the effect of twist averaging within the CE is to replace sums over the discrete set of single-particle 
orbitals by integrals over a volume of k-space. Consider, for example, a simulation cell of HEG containing an even 
number of electrons N. For each twist k s , the N/2 shortest Bloch vectors of the form k s + G s are doubly occupied. 
Integrating over twists therefore averages over the volume of k-space occupied by the first N/2 Brillouin zones of the 
simulation cell. The occupied region is a convex polyhedron that tends to the Fermi surface in the limit of infinite 
system size and has the correct volume at all system sizes. Since the sing le- particle KE k 2 /2 is a convex function of k, 
the small differences between the occupied region of k-space and the Fermi volume cause the CE twist-averaged HF 
KE to be slightly too large for finite systems. This systematic error, which exhibits visible shell-filling effects, decays 
with system size. 

Twist averaging within the grand canonical ensemble (GCE) also means taking the average of expectation values 
with respect to k s , but this time allowing the number of electrons to vary with k s . For any given k s , only those states 
that lie within the Fermi surface are occupied. This allows one to integrate over the Fermi volume in simulations 
with a finite number of particles, so that the HF KE of a HEG is exact at all system sizes. The KE at a given k s is 
obtained by summing the one-electron KE's of the occupied states. Values of k s with fewer occupied states therefore 
contribute less to the GCE average. 

We compare the efficiency of grid-based and Monte Carlo methods for integrating over the simulation-cell Bloch 
vector k s in Sec. Mil where we also discuss the use of the CE and GCE in HF calculations. 



E. Ewald interaction 



When simulating infinite periodic systems or finite systems subject to periodic boundary conditions, it is not 
possible to use the familiar 1 /r form of the Coulomb interaction because the sums over images of the simulation cell 
do not converge absolutely. The standard solution to this problem is to replace the Coulomb interaction by the Ewald 
interaction^ The 3D Ewald interaction is the periodic solution of Poisson's equation for a periodic array of point 
charges embedded in a uniform neutralizing background and is therefore appropriate for an electrically unpolarized, 
neutral system. Using the 3D Ewald interaction corresponds to adding a neutralizing background if necessary and 
calculating the Coulomb energy per simulation cell of a macroscopic array of identical copies of the simulation cell 
embedded in a perfect metal so that surface polarization charges are always screened* 2 * The Ewald energy for any 
particular electron configuration in a 3D system is 

where ve(j) is the Ewald interaction and vm = hm r _»o [ve(t) — 1/r] is the Madelung constant, which represents the 
interaction between a point charge and its own images and canceling background. These quantities may be evaluated 
efficiently using the Ewald formulas: 

1 v-^ 47rexp(- K 2 G 2 /2 + zG s • r) 2nn 2 ^ erfc [|r - R,|/(V2k)1 

VE ^ = nT, q2 ^ + E ^Tr^ ■ (6) 

1 4-7rexp[-G 2 /(4K 2 )] it t ^ erfc(K.R s ) 2k 



vm 



G s ^0 * R s #0 



where is the volume of the simulation cell. The value of the constant n does not affect Ve (r) or vm and may be 
chosen to maximize computational efficiency. The zero of potential has been chosen such that u_e(r) averages to zero 
over the simulation cell. The periodic function Ufi(r) has Fourier components 2 -^ ve{G s ) = Att/G 2 for G s ^ and 
ve(G s ) = for G s = 0. Setting k = l/(2y/e), where e is very small, Eq. Q gives 

1 ^ 47rexp(-eG 2 ) 1 1 47rexp(-eG 2 ) 1 f An exp(-e/c 2 ) 

"""n ^ g 2 = ^ ^ g 2 WF Loo P ' ( ) 



which will prove useful later on. 
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The analogous expression for the quasi-2D Ewald interaction is obtained by solving the 3D Poisson's equation for 
a 2D-periodic lattice of point charges subject to periodic boundary conditions in the plane and symmetric boundary 
conditions perpendicular to the plane, and is thus appropriate for planar and slab systems.— When evaluated in the 
plane of the charges, r± — 0, the 2D Fourier components VE(G s \\,r±) of the quasi-2D Ewald interaction VE(r\\,r±) 
are equal to 27r/G s || for G s m ^ and to for G s ii = 0. 



F. Structure factor and XC hole 



The analysis of the Coulomb and KE finite-size effects is most easily expressed in terms of the static structure factor 
(SF), the pair density, and the XC hole. The definitions of these quantities and relations between them are reviewed 
in this section. 

The SF is 

S(r, r') = | ([p(r) - p(r)}[p(r>) - p(r')}) = £ [</3(r)p(r')> - P {v)p{v')] , (9) 

where p(r) — ^ 5(r — r,) is the operator for the electron number density at position r, and p(r) = (p(r)) is its 
expectation value. In periodic systems, the Dirac delta functions are to be interpreted periodically: S[r— (r^ +R S )] = 
5(r — r,-). The SF is closely related to the pair density defined by 

P2 (r, r') = /j28(r r^(r' Tj )\ = ^S(r, r') + p(r)p(r>) S(r r')p(r'). (10) 

Another related quantity is the XC hole, p xc (r,r'), defined by 

p xc (r, r')p(r') = p 2 (r, r') - p(v)p(v>) = ^S(r, r') - S(r - r')p(r'). (11) 

Integrating Eq. pH|) with respect to r yields f Q p2(r,r')dr — (N — l)(X)j^( r ' ~ r j)) = {N — l)p( r ') and hence we 
obtain the sum rule J Q p X c(r, v')dr = — 1. The XC hole describes the suppression of the electron density at r caused 
by the presence of an electron at r'. 

It is often more convenient to work with the translationally averaged SF 

5(r) = l I 5(r' + r,r')dr', (12) 

and the analogous translationally averaged pair density pi (r) . These quantities have the periodicity of the simulation 
cell and may be expanded as Fourier series, the components of which are 

S(G S ) = ±{(p(G s )p*(G s ))-p(G s )p*(G s )}, (13) 
P2(G S ) = i (p(G s )p*(G s )) - ^ = ^ [S(G S ) 1] + Ip( Gs )p*(G s ), (14) 

where p(G a ) = J^i exp(— iG s • r^) is a Fourier component of the density operator^ Finally, the system-averaged XC 
hole is defined as 

Pxc(r) = i j p xc (r' + r, r')p(r') dr' = S(v) - S(v). (15) 
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G. Hartree and XC energies 

The Ewald interaction energy is the expectation value of the operator in Eq. 
Nv M 1 



Vew) = 




ll n j(H s ( r - r *w r ' - ^ ■ ) ) - r ') dv dv ' 



P2(r, r') ve{y — r') dr dr' (16) 
p xc (r,r>(r') [v E (r - r') - v M ] dv dv 1 + ~ / / p(r)p(r> B (r - r') dr dr' (17) 



n 



/n Jn ^ Jn Jn 



X - J2 MGs) [S(G S ) - 1] I + ^ E «b(G p )p(G p ) p *(G p ), 

1 G 3 #0 / G p #0 



(18) 



where use has been made of the sum rules J n p X c(r, r') dr = —1 and f n p(r') dr 1 — N. The first term in Eqs. (fTF)) and 
(TTH|) is the XC energy (the interaction of the electrons with their XC holes) ^ The second term is the Hartree energy 
(the interaction of the charge densities). The Hartree term in Eq. (|T5|) has been simplified by noting that p(r) has 
the periodicity of the primitive lattice and hence that p(G s ) vanishes unless G s £ {G p }. 

In practice the charge density and pair density converge rapidly with system size for interacting systems^ due to 
the fact that the XC hole falls off very quickly with r. For example, the nonoscillatory part of the XC hole falls off 
as r~ 8 for a 3D HEG. 1 - If the charge density is correct then the Hartree energy in a finite cell is exact, as can be 
seen from Eq. (|18[) : the Fourier components ve(G p ) are equal to 4ir/Gp and p(G p ) is proportional to the number 
of primitive cells in f2, so the Hartree energy per electron is independent of system size. The finite-size errors in 
the interaction energy given by Eq. I|17p must therefore be caused by the slow convergence of the Ewald interaction 
ve (r) — vm in the XC energy. 

A power expansion of the Ewald interaction about r = gives^ 



ve(t) -v m = - + —r^Wr + [-~), (19) 



1 2tt t / r 4 

~ + ^ r Wr + O — 

where the tensor W depends on the symmetry of the lattice. W is the identity matrix for a lattice of cubic symmetry. 
For large simulation cells the first term in the expansion dominates in the region where the XC hole is large, but for 
typical cell sizes the second term can be significant. Unlike the Hartree energy, we do not want the effect of periodic 
images in the XC energy: the interaction between each electron and its XC hole should just be 1/r. This is enforced 
in the MPC interaction^^ 

In HF theory, unlike QMC and reality the exchange hole is long-ranged (the nonoscillatory tail falls off as r~ 4 ) 
and the pair density is slowly convergent with system size^ This gives an additional source of finite-size error, even 
when the MPC interaction is used, as discussed in Appendix lAl 

H. MPC interaction 

The MPC interaction operator^i 4 - is 

Vmpc = \ Y1 f( Yl ~ r J') + H / { VE ( r * ~ r ) ~ f( ri ~ r )l dr 

~J^p(r)p(r'){v E (r-r')-f(r-r')} drdr', (20) 

where /(r) is 1/r treated within the minimum- image convention in the simulation cellj 1 ^ Assuming the pair density 
and the charge density have converged to their infinite-system forms, the MPC electron-electron interaction energy is 

(Vupc) = 1 [ [ p xc (r,r')p(r')f(r-r')drdr' + 1 f f p(v)p(r')v E (r - r') dr dr' . (21) 
x ' * Jn Jn z Jn Jn 

Hence the Hartree energy is calculated using the Ewald interaction while the XC energy is calculated using 1/r 
(within the minimum- image convention), as desired. The MPC interaction energy per electron therefore converges 
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more rapidly with system size than the Ewald interaction energy. One can avoid the need to know p exactly by 
replacing it with the approximate (but usually highly accurate) charge density p A from a DFT or HF calculation in 
Vmpc- The error due to this approximation is 0(p — pa) 2 - Comparing Eqs. (|2Tj) and (fTTj) . we see that the difference 
between the Ewald and MPC XC energies involves the operator (v E — v m — f), which vanishes as the size of the 
simulation cell goes to infinity. So the Ewald and MPC XC energies per particle are the same in the limit of large 
system size, even if an approximate charge density is used. In practice the first term of the MPC interaction is 
evaluated in real space, the second term is evaluated in k-space, and the third term is a constant; 19 

Vmpc - ^ E ^ ( r * ~ ^ + ^ E E MG P ) - /(G p )] p A (G p ) exp(zG p • r 4 ) 

+ f ~Nf oPAo - -L J2 IMG P ) - /(Gp)] P * A (G p )p A (G p ) + ^/opIoPao , (22) 

where /o and p A o are the G s = components of / and p A . Although /(r) has the periodicity of the simulation cell, 
its Fourier components are only required on primitive lattice vectors G p . These Fourier components are evaluated 
numerically in advance, a procedure that requires some care because /(r) diverges at r = and is nondifferentiable 
at the boundary of the Wigner-Seitz cell of the simulation cell. Once the Fourier components have been obtained, 
the MPC interaction is much quicker to evaluate than the Ewald interaction because (i) there is no real-space sum 
over lattice vectors and (ii) the k-space sum runs over primitive-cell G p vectors only, so the number of G p vectors to 
include in the sum does not grow with system size. 



I. Finite-size correction to the Ewald energy in 3D 



Assuming that the charge density (and hence Hartree energy) and the Fourier components of the SF converge rapidly 
with system size, it follows by comparing Eq. (TTS)) with its infinite system-size limit that the finite-size correction to 
the 3D Ewald interaction energy is 

AVew = Y 72^ / VE[k) [s(k) ~ 1] dk ~ si E Ve{Gs) [s{Gs) ~ 1] ~ VM ) ' (23) 

\ l ) Jk<oo G s ^O I 

where we have noted that vm — > as the system size tends to infinity. Since S(k) — > 1 as k — > oo, the sum and the 
integral converge, allowing us to include factors of exp(— ek 2 ) in the summand and integrand without affecting AVe w 
if e is small enough. Substituting for vm using Eq. 1(5)) then gives 




v E (k)S(k)exp(-ek 2 ) dk-~ E v E {G s )S(G s )exp(-eG 2 s )\ 

G 3 /0 / 



A^w « - | / v E (k)S(k) exp (-ek 2 ) dk - - > j v E (G s )S(G s ) exp (-eG 2 s ) \ . (24) 

'fc<oo 

The convergence factors are now required to keep the summation and integration finite, even though they do not 
affect the value of the expression as a whole. 

An obvious contribution to the finite-size correction is apparent from the form of Eq. (|24|) . In interacting electron 
systems with cubic (or higher) symmetry, S(k) = i]k 2 + 0(fc 4 ) for small k, where r\ is a constant.— The function 
v E (k)S(k) — 4TrS(k)/k 2 therefore tends to a well-defined limit as k — > 0, suggesting that much of the difference 
between the sum and the integral in Eq. (|24[) is caused by the omission of the G s = term from the summation. 
This argument leads to a finite-size correction of the form derived by Chiesa et al.<^ 

N AnSjk) 2nN V 

AVE ^mk™ a ^ 2 - = —- (25) 

Since (Ve w ) is proportional to system size, the relative error in the Ewald energy falls off as 0(iV _1 ). In a 3D HEG, 
the random phase approximation (RPA) implies that rj = l/(2w p ), where lu p — y/ 'iirN /Q = y3/Vf is the plasma 
frequenc y 20 ^ 1 giving^ 

AVew = ^. (26) 



7 



These approximate arguments may be made more precise and given an appealing physical interpretation as follows. 
According to Eq. (fl"5j) . S(r) = p xc (r) + S(r) can be viewed as the localized charge distribution of an electron at 
the origin and the system- averaged XC hole surrounding it. More precisely, because the simulation cell is repeated 
periodically, S(r) is a superposition of many such localized charge distributions, one centered in every copy of the 
simulation cell, i.e., 5*(r) = J2n S\ oc (r — R s )- If we assume that the XC hole is well localized within the simulation 
cell, which must be the case if S(k) has converged with respect to system size, this decomposition is unambiguous. 
It is then easy to show that 



S(G S ) = f 5 loc (r) ex P HG s • r) dr. 

J r<oo 



(27) 



The discrete Fourier components of the periodic function S(r) are therefore equal to the corresponding components of 
the continuous Fourier transform of the localized function S\ oc (r) . If S\ oc (r) is convolved with a very narrow normalized 
Gaussian (47re) -3 / 2 exp(— r 2 /4e) before the Fourier transform is taken, S(k) is multiplied by the convergence factor 

The convolution smears out the delta function slightly, but has no other discernible 

The integral is the value at the 



exp(— ek 2 ) appearing in Eq. ([24 
effect on the form of Si oc (r). 

We can now interpret the two terms between the large parentheses in Eq 
origin of the potential 



^loc,e(r) = / 

Jr'< 



sW( r ') 



dr' 



(28) 



corresponding to the aperiodic charge density Sio C ,e(r) obtained by convolving S\ oc (r) with the very narrow Gaussian. 
The summation [including the missing G s = term, which is well-defined for systems of cubic symmetry or if S(G S ) 
is replaced by its spherical average 2 ^] is the value at the origin of the potential 



4> e (r) = </>loc,e(r - R s ) 



R 



(29) 



of an infinite periodic lattice of copies of Si OCie (r). The sum rule f r<OQ Si oc ,e(r) dr = ensures that Si oc ,e( r ) nas 
no monopole and the system averaging of the symmetric function S(r, r') = S(r',r) ensures that Si oc ,e(r) has no 
dipolei 2 ^ If the system has cubic symmetry or we approximate S , i oc .e(r) by its spherical average as proposed in Ref. 0, 
the quadrupole vanishes too and <^i oc ,e(r) decays rapidly enough to ensure that the summation in Eq. (|29[) converges. 
Equation (|24|) can then be rewritten as 



N 



</M0) - 



5> loc (R s ) - i lim v E (k)S(k) 

R s 



N [ Air SCk) v-^ , .„ . . 



Q k^o k 2 



R s #o 



(30) 



a result that can also be obtained using the Poisson summation formula 24 (which we have, in effect, derived). The 
first term is the finite-size correction from Eq. (|25[) and the second term is small, as explained below. 

The nonoscillatory tail of the spherical XC hole of a 3D HEG is of the form p X c(^) = — Ar -8 , where A is a 
constant^ It arises from the 0(k 5 ) term in S(k)£2. The total XC charge lying further than r from the origin is 
therefore — 47rA/(5r 5 ), so 4>i oc (r) = 47rA/(5r 6 ) for large r. Hence 



N 
~2 



E 

R s #0 



N f°° 



4ttA 
5r f 



47r-r 2 dr = 0{N~ 1 ), 



(31) 



where Rq is the radius of a sphere of volume O. Thus, the remaining error in the XC energy per particle not accounted 
for by Eq. {25J is 0(N- 2 ). 

In inhomogeneous systems, /O xc (r) may not be spherical, causing 4>i oc (r) to decay more slowly at large r. In particular, 
if 5ioc(r) has a nonzero quadrupole moment, <fi\ oc (r) cx r~ 3 and the sum over R s fails to converge absolutely. The 
error not accounted for by the XC correction proposed by Chiesa et al^ is then of the same order as the correction 
itself. These additional errors are related to the behavior of S(]i)/k 2 near k = and are analyzed in Sec. IV Dl 

The MPC and XC correction methods are compared in Sec. IIVI The near equivalence of the MPC and the AVgw 
correction in cubic systems is proved very directly in Appendix [Bl 
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J. Finite-size correction to the KE in 3D 



According to the inhomogeneous generalization 26 ' 27 ! 28 of the Bohm-Pincs RPA, 29 which is believed to provide an 
accurate description of long-ranged correlations of electrons in solids, the wave function of a many-electron system 
may be approximated as 



* = *3 cxp £ u(G s )Ap*(G s )Ap(G s ) ] j 



(32) 



where Ap = p — p and ^ s has short-ranged correlations only. Expressed in terms of the coordinate operators, the 
RPA wave function takes the familiar 5 form 



* = * s cxp [ - u(ri -Tj)+Y, Xfa) j , (33) 

where x( r ) — — Jn u ( r ~ r ')p( r ') <^ r '- 

The long-ranged correlations are described by the function u(r), which has the periodicity of the simulation cell 
and inversion symmetry. At large r, u(r) is spin-independent and, in a 3D system, usually decays like 1/r. However, 
u(r) is necessarily restricted in a finite simulation cell, thereby biasing the KE. 

In a VMC simulation, the KE is evaluated as the average of the sampled values of 5 

f = --V 2 log(*) =f s -±-Y, U ( G ^ 2 [A/3*(G s )Ap(G s )] , (34) 

Gs^O 

where T s = —V 2 log(^ s )/4 and V = (Vi, .... Vat) is the 3A-dimensional gradient operator. It can easily be shown 
that V 2 [Ap*(G s )Ap(G s )} = -G 2 p*(G s )Ap(G s ) - G 2 p(G s )Ap*(G s ) + 2G 2 S N. Since (Ap(G s )) = 0, and hence 
(p{G s )Ap{G s )) = (Ap*(G s )Ap(G s )), it follows that 

T) = (Ts) + ^z E G 2 s [u(G s ) (Ap*(G s )Ap(G s )) - Nu(G s )] (35) 
f ") + §i E G 2 u(Gs)S*(G s ) E G 2 s u (Gs). (36) 

G 3 ^0 G 3 ^0 

We assume that (T s ) is exactly proportional to the system size (i.e., any finite-size error in (T s ) has been eliminated 
by twist averaging or the use of HF corrections) and concentrate here on the long-ranged finite-size errors arising from 
the Jastrow factor. Although the sum over G s in Eq. ([33]) converges, the two contributing terms in Eq. (f3"6")) diverge. 
As in the analysis of the Coulomb errors in Sec. Ill II this difficulty can be overcome by the inclusion of convergence 
factors, which are to be understood in the rest of this work. 

In practice u(k) has roughly the same form at different system sizes, since its small-k behavior is determined by the 
RPA.i Hence, in the infinite-system limit, the sum over G s in Eq. (|36|) can be replaced by an integral without changing 
the function u(k). For a symmetric system, u(r) = —A/r for large r, where A is a constant, 29 so u(k) = —AnA/lt 2 
at small k. Therefore lim^o fc 2 w(k) is finite, and the leading contribution to the finite-size error is the omission of 
the G s = term in the second summation in Eq. (|36p . The G s = term in the first summation is less important 
because 5(k) = 0{k 2 ). This argument leads to the finite-size correction proposed by Chiesa et al.: 1 

Nit A . . 

In the HEG, where the RPA implies that A = l/u> p 21 ' 29 this correction becomes AT = w p /4. 



III. COMPARISON OF TWIST-AVERAGING METHODS 



HF theory is the simplest framework in which twist-averaging methods can be compared. Very large simulation 
cells and twist samplings can be used, allowing the convergence with cell size and number of twists to be assessed 
reliably. Some of the finite-size errors that affect real interacting systems are not present in HF calculations, but twist 
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FIG. 1: (Color online) Convergence of the calculated KE per electron (top panel) and exchange energy per electron (bottom 
panel) of a 338-electron simulation cell of HEG at r s = 1 a.u. as a function of the number of twists for the three different CE 
twist-averaging methods described in the text. Because of the finite size of the simulation cell, the calculated KE and exchange 
energy do not converge to their exact infinite-system limits as the number of twists increases: the KE shows a small positive 
bias and the exchange energy a large negative bias. 



averaging is only intended to remove single-particle errors and the HF framework provides a valid test of how well it 
achieves this aim. 

The first issue is the choice of quadrature. The integrations over the simulation-cell Brillouin zone that yield twist- 
averaged energies cannot be carried out exactly and must be approximated by sums over finite sets of k s -points. We 
have considered three choices for the set of points: (i) a uniform Monkhorst-Pack grid 30 centered on the r-point of the 
simulation-cell Brillouin zone, (ii) a uniform grid centered on the Baldereschi-pointii of the simulation-cell Brillouin 
zone, and (iii) a random sampling within the simulation-cell Brillouin zone. All three choices yield identical results as 
the number of electrons N or the number of twists M tends to infinity, but the two limits are not equivalent: the fully 
twist-averaged (M — > oo) exchange energy depends strongly on N in both ensembles, while the fully twist-averaged 
kinetic energy depends weakly on N in the CE only. Since practical QMC simulations are unlikely to use very large 
simulation cells or numbers of twists (large numbers of twists are difficult because the full many-electron trial wave 
function must be constructed and stored for each twist), the rates of convergence with N and M are important. 

If the system is an insulator, the same bands are occupied at every k s and the integrand (e.g., the total kinetic 
energy as a function of k s ) is very smooth. The sampling theorem then ensures that estimates of the integral obtained 
using uniform twist grids converge very rapidly as the number of twists M is increased. If the twists are distributed 
randomly, the statistical error in the estimate of the integral decays more slowly, like M -1 / 2 . The most rapid 
convergence with number of twists and system size is obtained using a uniform grid of twists offset to the Baldereschi 
point 11 of the simulation-cell Brillouin zone. 

In metals, the integrand is discontinuous because of the sharp Fermi surface and the convergence with system size 
and number of twists is much slower. Figure [T] shows the HF kinetic and exchange energies of a face-centered cubic 
(FCC) simulation cell of HEG containing 338 electrons at r s = 1 a.u., calculated using sets of twists of various sizes 
generated in all three ways. As for insulators, energies calculated using random twist sampling converge slowly as the 
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FIG. 2: (Color online) Convergence of the calculated KE per electron (top panel) and exchange energy per electron (bottom 
panel) of a HEG at r s = 1 a.u. as a function of N for the three different CE twist- averaging methods described in the text, 
each with eight twists. The F-point and Baldereschi-point results have been offset for clarity. 



number of twists increases. The most rapid convergence is again obtained with a uniform Monkhorst-Pack grid of 
twists centered on the Baldereschi point of the simulation-cell Brillouin zone. The twists on a T-point Monkhorst-Pack 
belong to stars of symmetry-equivalent twists yielding identical energies. The symmetry can be used to reduce the 
number of trial wave functions that have to be constructed, optimized and stored per twist, but does not decrease the 
total number of Monte Carlo samples required to obtain a given statistical error and does not affect the conclusion 
that the Baldereschi-point grid is the most efficient. Because the simulation cell only contains 338 electrons, the KE 
and exchange energy do not converge to their infinite-system limits as the number of twists increases. The small 
positive error in the calculated KE is an artifact of the CE twist-averaging algorithm, as discussed in Sec. Ill Dt and 
disappears when GCE twist averaging is used. KE's in QMC simulations suffer from much larger finite-size errors 
due to long-ranged correlations (see Sec. Ill jp . but these are absent in HF theory. The large negative finite-size error 
in the exchange energy is not caused by the CE twist-averaging algorithm and is not removed by GCE averaging, but 
arises from the compression of the exchange hole into the simulation cell. 

Figure [2] shows the convergence with system size of the CE twist-averaged HF KE and exchange energies of a HEG 
at r s = 1 a.u. in a face-centered cubic (FCC) simulation cell, calculated using sets of twists generated in all three 
ways. To highlight the differences between the three methods, we have used only eight twists in each case. Energies 
calculated using the uniform grid of twists centered on T converge the most slowly because of the large fluctuations 
that occur as the size of the simulation cell increases and shells of symmetry-equivalent k s + G s vectors cross the Fermi 
surface. Energies calculated using a random sampling of twists converge more rapidly with system size (although less 
rapidly with number of twists). Yet again, the best approach uses a uniform grid of twists centered on the Baldereschi 
point of the simulation-cell Brillouin zone. 

Figure [3] shows the error in the twist-averaged HF KE calculated with a very large set of random twists, using both 
the CE and GCE. The systematic bias in the CE average disappears when GCE averaging is used, but the large 
fluctuations in the GCE results outweigh the bias for all but the smallest simulation cells. These fluctuations arise 
from the variations in electron number inherent in the GCE method. Most QMC simulations are likely to use many 
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FIG. 3: (Color online) System-size dependence of the twist-averaged KE per electron of a HEG at r 8 = 1 a.u. in the CE and 
GCE. All calculations used a random sampling of 5120 twists. 
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FIG. 4: Bias in the KE per electron of a HEG at r s = 1 a.u. as a function of N, calculated using 5120 randomly chosen twists 
in the CE. The power-law fit yields a bias proportional to N~ 1,32 . 



fewer twists, rendering the GCE fluctuations even worse, so CE averaging is the more promising method despite the 
bias. Figure 0] shows the bias in the CE-twist-averaged KE as a function of N. The power-law fit shows that the bias 
per electron decreases relatively slowly with system size, scaling roughly as N~ 4 / 3 , as noted by Lin et al& 

IV. COMPARISON OF THE MPC INTERACTION WITH THE FINITE-SIZE CORRECTION TO THE 

EWALD ENERGY 



If the XC hole can be assumed to have converged to its infinite-system form then both the MPC interaction and 
the finite-size correction of Eq. ([23]) are good solutions to the problem of finite-size effects in the XC energy of a 
cubic system. For low-symmetry systems the MPC interaction should continue to be a good solution, whereas the 
correction to the Ewald energy cannot be applied straightforwardly. On the other hand, if the simulation cell is too 
small to contain the infinite-system XC hole, but the SF is known analytically at small k, then this information can 
be included in the XC correction but not the MPC interaction, so the XC correction may work better. In practice the 
difference between the MPC energy and the corrected Ewald energy for cubic interacting systems is very small when 
the Ewald interaction is used to generate the configuration distribution, as demonstrated by the data shown for 3D 
HEGs at three different densities in Table U In each case the difference of MPC and Ewald energies is approximately 
equal to (but slightly greater than) AVew- 

It is shown in Appendix [5] that the long range of the exchange hole causes the MPC energy to be slowly convergent 
when the interactions are treated within the HF approximation. The finite-size correction constructed using the 
known small-fc behavior of the HF SF therefore performs better than the MPC interaction in HF calculations. 
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r s N (Emvc — -EEwaid)/^ (a.u. / elec.) AVew/A (a.u. / elec.) %age difference 



1 


54 


0.007 81(1) 


0.008 02 


2 6(1)% 


1 


102 


0.004 137(9) 


0.004 245 


2 5(2)% 


1 


226 


0.001 89(1) 


0.001 92 


1.6(5)% 


3 


54 


0.001 551(4) 


0.001 543 


0.5(3)% 


3 


102 


0.000 802(2) 


0.000 817 


1.8(2)% 


3 


226 


0.000 365(1) 


0.000 369 


1.1(3)% 


10 


54 


0.000 242(1) 


0.000 254 


4.7(4)% 


10 


102 


0.0001319(4) 


0.000 134 2 


1.7(3)% 


10 


226 


0.000 060 5(7) 


0.000 060 6 


0(1)% 



TABLE I: Difference between total energies per electron evaluated using the MPC and Ewald interactions in twist- averaged 
DMC calculations for 3D paramagnetic HEGs at three different densities. The Ewald energy is used in the branching factor 
in the DMC simulation, so that the configuration distribution appropriate for the Ewald interaction is used in all cases. The 
DMC time steps were 0.003, 0.03, and 0.3 a.u. at r s = 1, 3, and 10 a.u., respectively, and the target population was more than 
400 configurations in each case. Twist angles were sampled randomly. At each density it was verified that the DMC energy 
did not change when the time step was halved, the configuration population was doubled, and the number of post-twist-change 
equilibration steps was quadrupled. The finite-size correction to the Ewald energy [Eq. (|26|l ] is shown for comparison. 



A 



Ewald propagation MPC propagation 

-Bew /N (a.u. / elec.) Empc/N (a.u. / elec.) Ee™/N (a.u. / elec.) Empc/N (a.u. / elec.) 



54 -0.068 69(6) -0.06715(6) -0.06818(6) -0.067 36(6) 

102 -0.067 62(3) -0.066 82(3) -0.067 68(6) -0.06715(6) 
226 -0.067 06(4) -0.066 61(4) -0.066 85(5) -0.066 77(5) 

TABLE II: Total energies evaluated using Ewald and MPC interactions for 3D paramagnetic HEGs at r s — 3 a.u. The results 
were obtained in twist-averaged DMC calculations, as described in the caption to Table [I] The Ewald energy was used in the 
branching factor in the results labeled "Ewald propagation," while the MPC energy was used in the results labeled "MPC 
propagation" (i.e., the XC hole was appropriate for the Ewald and MPC interactions, respectively). -Bew and -Bmpc refer to 
the interaction (Ewald and MPC, respectively) used in the local energies that were averaged to obtain the DMC energy. 



By the variational principle, the expectation value of the MPC Hamiltonian with respect to the Ewald ground-state 
wave function is greater than the expectation value of the MPC Hamiltonian with respect to the MPC ground-state 
wave function. The MPC energy obtained using DMC with the Ewald energy in the branching factor is therefore 
likely to be overestimated, and vice versa. An example of this effect is shown in Table UT1 When the Ewald interaction 
is used in the branching factor, the difference between the MPC and Ewald energies is given by AVe w - However, when 
the MPC interaction is used, the difference is less than AVe w - These results suggest that the MPC interaction distorts 
the XC hole in a finite system, while the Ewald interaction gives a better shaped hole, although the interaction with 
the hole is not quite right. We have directly verified that this is the case for a HEG at r s = 3 a.u., as can be seen in 
Fig. \5\ The Ewald XC hole converges to its infinite-system form much more rapidly than the MPC hole. The likely 
reason for this behavior is that the MPC Hamiltonian does not include corrections for finite-size errors in the KE. 



V. NONANALYTIC BEHAVIOR AT k = 
A. Examples of nonanalytic behavior at k = 

The XC correction discussed in Sec. Ill II works well for interacting systems of cubic symmetry. In other cases, 
however, the theory cannot be applied straightforwardly. We give two examples. 

For a general interacting system, the SF at small k can be written as 5(k) = (l/2)k T IF'k for some tensor W . If 
the system has cubic symmetry then W is proportional to the identity matrix and hmk— »o S(k)/k 2 is well-defined. 
Otherwise, this limit is undefined and it is not possible to add the G s = term to the sum in Eq. (|24p. 

Within HF theory, S(k) = Afc + 0(k 3 ) at small k, where A is a constant^ The limit of S(k)/k 2 as k — > is 
therefore undefined. Again the approach discussed in Sec. Ill II cannot be applied. 
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FIG. 5: (Color online) System-averaged XC hole in a 3D paramagnetic HEG of density parameter r a = 3 a.u. at different 
system sizes relative to the XC hole in a 118-electron HEG with the Ewald interaction. The Ewald and MPC interactions 
were used to generate the configuration distributions. Twist-averag ed VMC and DMC XC holes p™°(r) and p™°(r) were 
calculated, and the final XC hole was obtained using the extrapolated estimate p xc (r) ~ 2p^ iC (r) — Pxc MC ( r )- The error in 
the extrapolated estimate is second order in the error in the trial wave function.— 



B. Removing the problematic part of the SF 

Suppose that S(k)/k 2 is singular or otherwise ill-defined at k = 0, but that its smali-k behavior is known and is 
roughly independent of N. We can then introduce a model "structure factor" <S(,(k) that incorporates the nonanalytic 
behavior and define 5 a (k) = S'(k) — 5f,(k), so that linik^o S'a(k)/fc 2 is well defined. Starting from Eq. ((24)) and 
applying the Poisson summation formula^! to terms involving S a only yields 

(38) 

where 5 tt (r) is a localized charge distribution analogous to Si oc (r) and all convergence factors have been omitted. 
Since the k — > behavior of S a Qs)/k is known, and provided that Sf,(k) has a simple enough form, all three terms 
within the large parentheses in Eq. ([55)) can be evaluated straightforwardly. Moreover, since 5 a (k) is well-behaved as 
k — * 0, 5 (r) lacks the long-ranged tail present in S\ oc (r); the summation in the final term on the right-hand side of 
Eq. ()38|) therefore converges rapidly and should be small. This term is omitted from the approximate expressions for 
the finite-size correction obtained below, and therefore represents the error in these approximations. 

The finite-size correction obtained by evaluating all except the final term on the right-hand side of Eq. (|3"5)) is 
accurate when S^k) = S'(k) — S^k) is smooth, implying that S a (r) is short ranged. The model structure factor 
S'b(k) should therefore match the nonanalytic behavior of 5(k) as closely as possible. It is also sensible, although 
less important, to ensure that S(k) — S'b(k) is small. In practice, although S(k) — > 1 as k — > oo, the correction is 
most easily evaluated if Sb (k) — > as k — -> oo. A natural way of accomplishing this is to include a Gaussian function 
exp(— ak 2 ) as a factor. The parameter a should be small enough that the Gaussian changes little on the scale of 
the Fermi wave vector. In fact, although the reciprocal space summation and integration diverge in the a — * limit, 
their difference converges rapidly. One can therefore maximize the smoothness of 5a (k) by decreasing a until the 
calculated value of the correction has converged. 

A plausible alternative method- for dealing with leading-order nonanalyticities in S(k)/k 2 at k = is to replace the 
missing G s = term in the sum over G s in Eq. (fT5)) with an integral of ^(fc)5(k) over a sphere of volume (27r) 3 /fi. 
This approach may be cast into the framework discussed above by choosing «S&(k) = S(k)Q(Q — k), where Q is the 
radius of the sphere of volume (27r) 3 /f2 and Q(Q — k) is a Heaviside step function. The function 5 a (k) — S(k) — S*b(k) 
is then zero at the origin, so the first term inside the parentheses in Eq. (|3"5)) vanishes. Unless the lattice is very 
asymmetric, St(G s ) is zero for all nonzero G s , and the third term inside the parentheses in Eq. (I38p also vanishes. 
Hence 
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FIG. 6: (Color online) HF Ewald exchange energy per electron of a 3D paramagnetic HEG of density parameter r s = 1 a.u. 
against particle number N with k s = ("PBC") and twist- averaged boundary conditions within the CE ("TABC"). The 
corrections of Eqs. (| A5|) and (|41|) have been applied to the data labeled "TABC (corr. A)" and "TABC (corr. B)," respectively. 
An FCC simulation cell is used. The uncertainty in the twist-averaged data due to the use of a finite number of twist angles is 
small compared with the difference between the twist-averaged data and the infinite-system result. 



In this case, however, the sharp cutoff in S&(k) leads to slowly decaying oscillations in Sb(r) and therefore 5 a (r). 
These oscillations fall off as r~ 2 and can never be regarded as negligible. Unless S(k)/k 2 is constant for k < Q, in 
which case this correction is accurate by construction, the neglected real-space term in Eq. (f3"TJ)) is of the same order 
as the correction itself. 



C. Finite-size corrections in HF theory 

Suppose S(k) = Xk + 0(k 3 ), as is the case for systems of cubic symmetry in HF theory. The divergence of S(k)/k 2 
as k — > prevents Eqs. (j2"5|) and ([50]) from being used to obtain finite-size corrections. Let S&(k) = Afcexp(— ak 2 ). 
Working in the a — > limit, Eq. ([55]) becomes 

\ G s ^O / 

where C HF = 2.8884, 2.8372, and 2.8882 for FCC, simple cubic (SC), and body-centered cubic (BCC) simulation cells, 
respectively^ and we have noted that the 0(k 3 ) term in 5 a (k) = 5(k) — Sb(k) causes S a (r) to fall off as r -6 , giving 
the 0(JV -1 / 3 ) correction. For a 3D paramagnetic HEG^ 7 A = (3/4)[fi/(3A r 7r 2 )] 1 / 3 , so 

AVfa-^^^ + OC^). (41) 

An alternative real-space treatment of HF finite-size errors can be found in Appendix [X] As shown in Fig. [6l 
both the real- and reciprocal-space approaches account for most of the HF Coulomb finite-size error, although the 
reciprocal-space approach performs better because it completely removes the 0(N 1 ' 3 ) error. 

D. Finite-size errors in the XC energy of low-symmetry systems 

For a general interacting system the SF can be written as 

i 

S ( k )=Y, J2 S lm(k)k l Y lm (6 k ,^), (42) 

even / m — — l 

where 9^ and (f>^ are the polar and azimuthal angles of k and Yi m is the (I, m)-th spherical harmonic. The odd-l 
components are zero by inversion symmetry. Guided by the RPA, we assume that S(k) is quadratic near k = 0, 
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and hence that Soo{k) oc k 2 . If the quadratic form is nonspherical, however, I = 2 components are also present and 
linik^o S(k)/k 2 depends on the direction in which the limit is taken; there is then a point discontinuity at k = 0. 
Equivalently, the I = 2 component gives rise to the quadrupole moment in Si oc (r), which leads to the additional errors 
discussed in Sec. Ill II 
Let 

2 

S b (k)= S 2m (O)fc 2 Y 2m (0 k ,0 k )ex P (-afc 2 ), (43) 

m=-2 

and S a (k) = S(k) — S&(k), where a is such that Sb (k) is long-ranged in k-space compared with the Fermi wave vector. 
Applying Eq. (|38|) and taking the limit a — > 0, we find that 

A^ew = f I £ ^ 2m (0) £ y 2m (^> G j) +0(iV-V3). (44) 

y m=- 2 g s5 ^o y 

In particular, it can be seen that the O(N ) finite-size correction obtained using the spherically averaged SF is 
incomplete, and that there is in general another correction of O(N ) due to the low symmetry of the simulation 
cell and the existence of the I = 2 component. If the XC hole has spherical symmetry, the extra correction is zero 
regardless of the shape of the simulation cell; if the XC hole does not have spherical symmetry, but the simulation 
cell does have cubic symmetry, the extra correction is again zero. Hence, if one is simulating a low-symmetry system, 
it is advisable to choose a simulation cell that is as close to cubic as possible. If this is not possible then one could 
evaluate the I = and I = 2 components of 5(k) at k = 0, and use Eq. (j4~4"| to compute the correction. The 0(7V -1 / 3 ) 
error in Eq. (|44[) arises from an assumed nonanalytic 0(k 3 ) term in S a (r). 



VI. HIGHER-ORDER CORRECTIONS TO THE KE 
A. Need to include higher-order corrections 

The need to include higher-order finite-size corrections to the KE is demonstrated in Fig. which shows the size 
dependence of the DMC energy of the 3D HEG. The XC- and KE-corrected Ewald data and the KE-corrected MPC 
data are in good agreement with each other, as expected. At low density (r s = 10 a.u.) the corrected data are 
almost independent of system size, indicating that the finite-size correction formulas are working well. However, at 
intermediate (r s = 3 a.u.) and high density (r s — 1 a.u.) it is clear that the QMC data are overcorrected when only the 
leading-order KE correction is applied. Since the finite-size correction AVe w to the interaction energy has been shown 
to be accurate, the problem must lie in the KE. It is clearly necessary to go beyond leading order when correcting 
the KE at intermediate and high densities. 

The Poisson summation formula can be used to demonstrate that higher-order terms are more important in the 
KE than the Ewald energy. If we assume that the XC hole is well localized within the simulation cell and that 
linik^o k 2 u(k) exists, the finite-size correction to the KE may be obtained from Eq. (j36)) as 

AT = ^ ( — / /,^k)[^k)-i],/k-- >_ r;;;,(G..[.S'(G.,!-:i; j (45) 



y r s #o / 



(46) 



where we have used the Poisson summation formula 24 and 



L(r) = -V 2 j u{r - r> xc (r') dr' (47) 

J r'<oo 

is the inverse Fourier transform of fc 2 w(k)[S'(k) — 1]. 

The leading-order behavior of the two-body Jastrow factor of a HEG at small k is^ 2 - u(k) = —Air (A/k 2 + B/k) 
within the RPA. Hence, at large r, u(r) = —A/r ~ 2B/(nr 2 ) and so V 2 u(r) = — 4S/(7rr 4 ) for r ^ 0. The finite-size 
correction to the KE is therefore 

_ NnA NB v f Pxc (r)dv 

AT - — -—t&L- pt^f- (48) 
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FIG. 7: (Color online) DMC total energy per electron against the reciprocal of the number of electrons in 3D paramagnetic 
HEGs of density parameter r s = 1 a.u. (top panel), r a — 3 a.u. (middle panel), and r s — 10 a.u. (bottom panel). The 
simulation parameters were as described in the caption of Table [T] The Ewald energy per electron is corrected by the addition 
of (AVew + AT A )/N = lo p /(2N) ["Ewald (corr.)"], while the MPC energy is corrected by the addition of AT A /N = u p /{4N) 
["MPC (corr.)"]. The corrected Ewald and MPC energies are hard to distinguish because they lie almost on top of each other. 
The higher-order KE corrections described in Sec. I VI Bl ( ATb /N plus the single-particle correction) are included in the data 
sets labeled "(new corr.)". 
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The first term is the correction of Eq. (|37p . while the second term gives an additional correction that falls off slowly 
as TV" 1 / 3 . So, even in the case of the HEG, where the next-to-leading-order correction to the Ewald energy falls off 
as -/V -1 , higher-order corrections to the KE may be important. 

The additional KE correction is due to the discontinuous gradient of k 2 u(k) at k = 0. A similar approach to that 
developed in Sec. IV Bl can be used to eliminate the leading-order nonanalytic contributions to the long-ranged part of 
V 2 u(r). Define F(k) = fc 2 u(k)[S(k) - 1] and write F(k) = F a (k) + F b (k), where F b (k) contains the O(k) contribution 
to — k 2 u{k) [as well as any anisotropic OQc ) terms], and is smooth and long-ranged in k-space. Then 



N I 1 
AT = — I — limF„ k 



fc<oo 



F b (k) dk 



G a ^0 



F b {G s ) 



N 
1 



E ^( R « 



(49) 



Note that, as shown in Sec. Mil the bias due to residual CE twist-averaged single-particle KE errors also falls off as 
TV^ 1 / 3 . If we include higher-order corrections for the neglect of long-ranged correlations, we should also correct for 
the residual error in the twist-averaged energy. 



B. Higher-order KE corrections 
Gaskell 3 - 2 - has derived the following expression for the small-fe limit of u(k) for the 3D HEG within the RPA: 



u(k) 



where 




f v E {k)N Y 
AS 2 (k) 1 \ Slu p J 
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1 /3 

is the HF SF, /cf ct = (67r 2 A r cr /J7) is the Fermi wave vector for particles of spin tr, N a is the number of particles of 
spin <7, and 
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(53) 
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where £ = (jVj — N\)/N is the spin polarization. 

Let F b (k) = 4nBk exp(— ak 2 ). This satisfies the requirements for F b given in Sec. IVI Al provided a is small. Then, 
by Eq. (g^) in the limit a -> 0, 
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N / Air A 
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4ttB 
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2/3 



Q{N- 2/3 ) = AT A 



AT F 



(55) 

(9(A^- 2/3 ), (56) 



where C 3D = 5.083, 5.264, and 5.086 for FCC, SC, and BCC simulation cells, respectively. The 0(A^" 2 / 3 ) error arises 
from the 0(r -3 ) term in u{r) at large r. 

The relative importance of the corrections for typical system sizes at three different densities is shown in Table [ 
The residual CE twist-averaged single-particle KE error is generally greater than ATg. This error can be estimated 
within HF theory! 3 - 3 - For real systems, the "infinite-system" HF energy would have to be evaluated in a large, finite 
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KE correction (a.u. / elec.) 



r s (a.u.) 


N 


SP corr. 


ATa/N 


AT B /N 


1 


54 


-0.002 8 


0.008 


-0.0016 


1 


130 


-0.000 65 


0.003 33 


-0.000 48 


3 


54 


-0.000 31 


0.00154 


-0.00017 


3 


130 


-0.000 072 


0.000 641 


-0.000 054 


10 


54 


-0.000 027 


0.000 254 


-0.000 015 


10 


130 


-0.000 006 


0.000 105 


-0.000 005 



TABLE III: Magnitude of different components of the finite-size correction to the KE per electron of a 3D paramagnetic HEG 
at different density parameters r 3 and system sizes TV in an FCC cell. The correction for residual single-particle errors after 
twist averaging in the CE ("SP corr.") is estimated as the difference between the infinite-system HF KE and the twist- averaged 
HF KE for the finite system. 

calculation. The effect of adding higher-order corrections (including the correction for the residual single-particle 
error) to the energy of a 3D HEG is demonstrated in Fig. \7\ The finite-size behavior of the QMC data is clearly 
greatly improved at r s = 1 and 3 a.u. 

For real systems we do not usually have an analytic result for the small-fc behavior of u(k). However, we have 
flexible forms of u(r) that can be optimized within QMC. By fitting a suitable functional form to the QMC-optimized 
u, we can extrapolate to the k = limit. We suggest that Eq. (|5ip be fitted to the QMC u(G s ) at the first two stars 
of nonzero G s vectors, and that Eq. (|55| should then be used to evaluate the KE correction. 

C. Low-k behavior of the Fourier-transformed two-body Jastrow factor 

The Fourier transform of the two-body Jastrow factor of a 3D paramagnetic HEG at r s = 3 a.u. is shown in Fig. 
[H The Jastrow factor consisted of polynomial and plane-wave expansions in electron-electron separation, 3 4 which 
were optimized by variance minimizatio n 3 ^ 36 followed by energy minimization^ As expected, the form of u(k) is 
largely independent of the number of electrons, and the small-fc behavior is well-described both by the RPA of Eq. 
([50)) and by the first two terms of the power-series expansion, Eq. (|51"j) . The RPA expression for the two-body Jastrow 
factor does not satisfy the Kato cusp conditions 3 ^ and hence becomes unreliable at large k. The small-A; behavior of 
u(k) for 54- and 226-electron HEGs is shown in Fig. [5] ft can be seen that the fit of the two-term RPA expansion 
to the QMC-optimized Jastrow factor is a fairly good approximation to the analytic RPA form within a sphere of 
volume (27r) 3 /f2, but that the fitted three-term RPA expansion is badly behaved, because one is simply fitting to the 
noise in the u(G s ) data. This is reflected in the corresponding results for the KE correction shown in Table [TV] The 
corrections obtained with the fitted two-term expansion are close to the analytic KE correction (leading-order and 
next-to-leading order terms). The leading-order correction can be thought of as being calculated on the assumption 
that k 2 u(k) is constant over the integration regions shown in Fig. [5J which is clearly inappropriate, and leads to the 
overcorrection for N = 54 electrons. The fitted three-term RPA expansion also gives an overcorrection. For HEGs, 
the analytic results given in Sec. IVIBl should of course be used. 

VII. FINITE-SIZE CORRECTIONS IN 2D SYSTEMS 
A. XC energy in 2D 

Consider a 2D-periodic system with simulation-cell area P. For a sufficiently symmetric system, S(k) = 7/c 3 / 2 + 
0{k 2 )M- Hence lim k ^ v E {k)S(k) = 27rlim k ^ S(k)/k = 0. So the 2D analog of Eq. (30]) is 

AV Evr = ~ ^oc(Ra). (57) 

R s #0 

In a 2D HEG the nonoscillatory XC hole is relatively long-ranged due to the reduced screening, decaying as p xc (r) = 
— Ar~ 7 / 2 , where A is a constant.— Hence the XC charge outside radius r is — 4-7rA/(3r 3 / 2 ) and the leading (monopolar) 
contribution to 0i oc (r) is proportional to r~ 5 / 2 at large r. [The dipole moment of the electron and its XC hole is zero, 
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FIG. 8: (Color online) Fourier transform of the fully optimized two-body Jastrow factor for a paramagnetic 3D HEG at r s — 3 
a.u. and different system sizes. For comparison, the RPA expression of Eq. (|5Up and the two-term RPA expansion [Eq. (I51|l ] 
are shown. 



Method N A (a.u.) AT (a.u.) 



Analytic RPA 


54 


3 


0.073 9 


Analytic 1-term exp. 


Any 


3 


0.083 2 


Fitted 2-term exp. 


54 


3.22 


0.079 1 


Fitted 3-term exp. 


54 


4.41 


0.097 6 


Analytic RPA 


226 


3 


0.077 5 


Fitted 2-term exp. 


226 


3.15 


0.0810 


Fitted 3-term exp. 


226 


3.53 


0.086 9 



TABLE IV: Finite-size correction to the KE of a paramagnetic HEG of density parameter r s — 3 a.u., calculated using different 
two-body Jastrow factors. The value of A = 3 corresponds to l/u p ; see Eq. (|54[) . We compare analytic results with those 
obtained by fitting to the QMC-optimized two-body Jastrow factors shown in Fig. [9] The "Analytic RPA" form is that of Eq. 
(|50|) . the "3-term exp." is that of Eq. (15 If) with an extra term — 4nC, the "2-term exp." is Eq. (|51|l . and the "1-term exp." is 
Eq. J51J with B = 0. 



while the quadrupole^ contribution to 0i oc (r) is proportional to r 3 .] Hence 

AVW oc ~ Yl R * 5/2 = OiN- 1 ' 4 ), (58) 

R s #0 

since the length of every simulation-cell lattice vector R s appearing in the summation is proportional to yN. Unlike 
the 3D case, therefore, AVew — * as N — ► oo. This conclusion was also reached, using a different approach, by Wood 
et alM 

To obtain the leading-order correction to the XC energy, we use the method of Sec. IV Bl Let Sb(k) — 
-yk 3 / 2 exp(-ak 2 ). Then, by the 2D analog of Eq. f38]). 



N 

2 1" ' (2n) 2 j — " P 



AFew = — I °+ 77TZTI / v E {k)S b {k) x 2irkdk 



p v E (G s )S b (G s ) \ +0(N- 1 / 2 ) 



G s ^0 / 

= ^+0(N-V>), (59) 

where C2D = 3.9852 and 3.9590 for square and hexagonal cells, respectively, and the a — > limit was taken in the final 
step. The 0{N~ 1 / 2 ) error is due to the quadrupole moment of S a (j) — S(r) — 5&(r). For a 2D HEG, 39 7 = 2~ 3 / i r^ 1 ^ 2 . 
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FIG. 9: (Color online) Two-body Jastrow factor for a paramagnetic 3D HEG at r s — 3 a.u. with 54 electrons (top panel) 
and 226 electrons (bottom panel). The Jastrow factors are as in Fig. [5] The RPA expression of Eq. (|50p is plotted, as are 
fits of the two- and three-term RPA expansions [Eq. (|51|l and Eq. (|51|) with an extra term — 47rC]. The fits were made to the 
QMC-optimized u at the first two and first three stars of G s vectors, respectively. The dotted line indicates the radius of the 
sphere whose volume is (27r) 3 /f2. 



Hence 
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This correction falls off very rapidly with r s 



2nr^ \ nN 



B. KE in 2D 



(60) 



For a symmetric 2D-periodic system^ u(k) — —ak 3 / 2 + 0(k - 1 ) and S(k) = 7A: 3 / 2 + 0(k 2 ). Hence, proceeding 
as in Sec. EDS we have F(k) = k 2 u{k)[S{k) - 1] = ak 1 / 2 + 0{k). Let F b {k) = ak 1 ' 2 exp(-afc 2 ). Then, by the 2D 
analog of Eq. (|4^|) . 
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where the a — > limit was taken in the final step. 
For a 2D HEG the HF SF isH 
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FIG. 10: (Color online) Twist- averaged DMC energy per electron against system size for a 2D paramagnetic HEG of density 
parameter r s = 20 a.u. The trial wave function was of Slater- Jastrow-backflow form,— the target population was 1536 
configurations, and the DMC energies have been extrapolated to zero time step. In each case the Ewald interaction was used 
in the DMC branching factor. "Ewald" and "MPC" indicate the interaction used to calculate the local energies. AT is given 
in Eq. (|64[) . respectively, and "SPC" denotes the single-particle correction to the KE (the difference of the infinite-system and 
CE twist-averaged finite-system HF KE's). 



where the Fermi wave vector for electrons of spin a is kp a = \J AirN a / P. The small-fc limit of the two-body Jastrow 
factor within the RPA is^i 
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(64) 



For real systems, we suggest that the Fourier transform of the two-body Jastrow factor be fitted to u(k) = —a/k 3 ^ 2 — 
b/k using the first two nonzero stars of simulation-cell G s vectors. Equation (ffTTj) should then be used to calculate 
the KE correction. 



C. Effectiveness of 2D KE correction 



We illustrate the effectiveness of the KE corrections in a 2D HEG at low density in Fig. [TOl The XC correction 
[Eq. (f6"D|) ] is negligibly small at this density. However it is clear that applying finite-size corrections to the KE alone 
is not sufficient to obtain accurate results. The MPC interaction gives significantly smaller finite-size errors than the 
Ewald interaction; nevertheless it is clear that extrapolation is necessary. 



VIII. FORMULAS FOR FINITE-SIZE EXTRAPOLATION 



A. Finite-size extrapolation 



In nearly all QMC studies of condensed matter to date it has been necessary to extrapolate energy data to infinite 
system size by means of an assumed relationship between energy and particle number. These formulas contain free 
parameters, including the infinite-system energy, which are determined by a fit to the QMC data. Despite the existence 
of sophisticated methods for treating finite-size errors, it is likely that some form of extrapolation will continue to be 
necessary for accurate work. In this section we analyze the performance of fitting formulas that have been proposed 
in the literature and consider how best to extrapolate QMC energies to infinite system size. 

Throughout this section we denote the QMC energy per electron of an ./V-electron system by e(N) and we denote 
the HF energy, KE, and interaction energy per electron by enp(N), tnp(N), and vn F (N), respectively. We assume 
that the same k s is used in both the QMC and HF calculations (or that twist averaging is applied in both cases). 
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B. Finite-size extrapolation formulas for the HEG 



The exact size-dependence of the HF energy of the fluid phase of the HEG is 

e HF ((») = ejjF (N) + At HF (N) + Av KF (N) , (65) 

where AtuF(N) = £hf(oo) — £hf (-AT) and Avhf(N) = i>hf(oo) — vhf(N). The forms of At^F(N) and Avrf(N) for 
a 3D paramagnetic HEG can be seen in Fig. [T] Both are oscillatory functions of N due to single-particle finite-size 
errors. The fluctuations in the exchange energy and the KE are strongly correlated, although those in the KE are 
larger. For further discussion of the single-particle finite-size errors in HF theory, see Sec. Ill Dl and Ref. 6. There is 
also a systematic error in the HF exchange energy, caused by the compression of the exchange hole, as discussed in 
Appendix [XI For a Wigner crystal, Ceperley^, suggested the fitting form 

e (oc) = e (Ar) + - f7 ^_ ) (66) 

where d is the dimensionality and c is roughly independent of r s . This is consistent with the form of the 3D XC 
correction [Eq. (j25|) ] and the leading-order correction to the KE [Eq. ([37]) ]. For an interacting Fermi fluid, CeperleyS 
suggested that the HF extrapolation is appropriate at small r s , while the Wigner-crystal extrapolation is more 
reasonable at large r s . He therefore proposed using an interpolation of Eqs. ((65)) and (|66|) . 



/ 1 JV -3/d 7 .3/2 \ 

e(oo) = e( J yO + A to (JV) + ^ s ^ + — . (67) 

For their study of the 3D HEG, Ceperley and Alder 44 used the two-parameter form 

e (oo) = e(A0 + aAi H F(A0 + ^7d, (68) 

where a and c are fitting parameters that vary with density. The parameter a may be thought of as the ratio of the 
actual electron mass to the effective mass within Fermi liquid theory. One therefore expects a « 1 in weakly correlated 
systems. Alternatively one can estimate a = At(N) / Atnp(N) . The parameter c accounts for the Coulomb finite-size 
effects in the XC energy and the neglect of long-ranged correlations in the KE. This form has also been used for the 
2D HEG, 41 although our analysis (see Sec. IVIIj) suggests that a term of the form cN~ 5 / 4 would be more appropriate 
than cN- 3 / 2 . In their studies of the 3D HEG, Ortiz et al& tested both Eqs. JBSJ and (JBTJ). They found that the two 
formulas give very similar results, but that c in Eq. (fBT)) was a strong function of r s . 



C. Comparison of extrapolation formulas 

Consider the extrapolation formula 

e(oo) = e(N) + aAt HF (N) + bAv HF {N) + ^- (69) 

for a 3D system, where a, 6, c, and 7 are parameters to be determined by fitting, which are allowed to vary with 
density. Imposing the constraint 6 = and 7=1 gives Eq. (f6"8"|) . The results of fitting Eqs. I|69p and (f6"T)l to DMC 
data for paramagnetic Fermi fluids at r s = 1, 3, and 8 a.u. are shown in Table fVl 

The extrapolated energies can be compared with the infinite-system limit of the Slater- Jastrow DMC energies 
obtained using twist averaging at r s = 1 and 3 a.u., as shown in Fig. [7] and quoted in the caption to Table IVl In each 
case the optimal value of b is approximately 0, and the x 2 value does not increase greatly when b = is imposed. 
Setting a = b (i.e., using the HF total energy to extrapolate away single-particle finite-size errors) gives a very poor fit 
to the data and introduces significant bias into the extrapolated energy. Both of these effects are caused by the slowly 
decaying systematic error in v-hf(N) due to the long-ranged exchange hole; this error does not have a counterpart in 
the QMC data to which the formula is fitted. At high densities the fit can be improved considerably by allowing 7 to 
vary; however the extrapolated energies are then biased. It is preferable to impose the known behavior 7 = 1. Setting 
the effective mass a equal to 1, which is also implicit in Eq. (ftJT)) . greatly increases the x 2 value of the fit, but does 
not significantly bias the extrapolated energy, because it simply reduces the amplitude of the oscillations in the fitted 
energy. Using Eq. (|67p or Eq. (|69[) with a = 1 is unreliable with small numbers of data points, however. Furthermore, 
Eq. (|67p is likely to be poor at low density because of the inclusion of Avhf^). Note that where the fits are good 
the effective mass ratios a are in good agreement with one another, and they increase with r s . 
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TABLE V: Results of fitting Eqs. (|69jl and (|67jl to DMC energy data for a paramagnetic Fermi fluid. Seven different system 
sizes (N = 18, 54, 118, 226, 338, 458, and 566) were used for each density, and the statistical error bars in the total energy were 
around 0.00001-0.00005 a.u. per electron. Time steps of 0.0015, 0.033, and 0.1 a.u. were used in the simulations at r s = 1, 3, 
and 8 a.u., respectively. The target population was 3200 configurations in each case. The wave-function was of Slater-Jastrow 
form, i.e., backflow was not used. The constraint a = 6 leads to an enormous \ 2 value at r s = 8 a.u. The infinite-system 
DMC energies from the twist-averaged DMC calculations at r s = 1 and 3 a.u. are 0.5880(6) and —0.06623(3) a.u. per electron, 
respectively. 



In summary, if single-particle finite-size errors are to be removed by extrapolation using Eq. (|69|) then only the HF 
KE should be used in the extrapolation formula (i.e., b should be 0), and some attempt should be made to compute 
the effective mass a. In 3D the exponent 7 should be 1, while in 2D it should be 5/4. 

IX. SIZE-DEPENDENCE OF BIASES IN DMC ENERGIES 

Figure QT] shows that the time-step bias in the DMC energy per particle has nearly the same form over the range of 
system sizes typically encountered in DMC simulations. A time step judged to be accurate in a small system should 
therefore continue to be accurate in a larger system. To exaggerate the bias, most of the results shown in Fig. fTTI were 
obtained using a simple Slater trial function with no Jastrow factor; the bias is greatly reduced if a more accurate 
trial wave function is used, as can also be seen in Fig. 111! 

For any given system, the DMC population-control bias should fall off roughly as Nq , where Nc is the target 
population, 46 so we have plotted the DMC energy against N^ 1 in Fig. |T2j Unlike time-step bias, population-control 
bias grows with system size. However, the increase in the bias with system size is slow. Population-control bias is 
caused by the correlation of fluctuations in the local energy and the DMC branching factor J 6 . Fluctuations in the local 
energy increase as N 1 / 2 . If the exponential branching factors can be approximated by the first two terms in the Taylor 
expansion of the exponential then fluctuations in the branching factor increase as N 1 ^ 2 . So the population-control 
bias in the energy per particle is roughly independent of system size. However, the fluctuations in the exponential 
branching factor grow more rapidly than N 1 / 2 in large systems, causing the bias to increase. Improving the accuracy 
of the trial wave function reduces population-control bias, as can be seen in the upper panel of Fig. 1121 

X. CONCLUSIONS 

We have carried out a detailed study of finite-size effects in QMC calculations and have described a number of 
approaches for reducing or correcting them. Twist averaging greatly reduces the magnitude of single-particle finite- 
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FIG. 11: (Color online) DMC energy per electron against time step for a paramagnetic 3D HEG ol density parameter r 3 — 4 a.u. 
at various system sizes N. A Slater wave lunction was used, except for the one curve labeled "Jas," in which a Slater-Jastrow 
wave lunction was used. Twist averaging was not applied. 
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FIG. 12: (Color online) Upper panel: DMC energy per electron against reciprocal of target population for a paramagnetic 
3D HEG of density parameter r s — 4 a.u. at various system sizes N. Lower panel: gradient of the population-control bias 
(derivative of the DMC energy per electron with respect to the reciprocal of the target population) against system size. The 
DMC time step was 0.03 a.u., and a Slater wave function was used, except for the one curve in the top panel labeled "Jas," in 
which a Slater-Jastrow wave function was used. Twist averaging was not applied. 
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size errors, although residual single-particle errors due to the wrong shape of the CE twist-averaged Fermi surface are 
still significant in studies of the HEG. One can calculate these errors within HF theory, and hence correct for them. 

Finite-size effects in the XC energy should be eliminated, either by adding a correction to the Ewald energy or by 
using the MPC interaction to calculate the final energies (although the Ewald interaction should be used to generate 
the configuration distribution, since the MPC interaction distorts the XC hole in finite systems). Finite-size corrections 
must also be applied to the KE. For HEGs, where analytic expressions for the \ow-k behavior of the two-body Jastrow 
factor are available, we have found that it is important to include both the leading- and next-to-leading-order KE 
corrections at intermediate and high densities. The resulting QMC energy data are almost independent of particle 
number at typical system sizes. For real systems we recommend fitting the QMC-optimized Jastrow factor to Eq. 
(f5Tj) at small k, then using Eq. ([55| to compute the correction to the KE. 

Within HF theory the long-ranged nature of the exchange hole leads to additional errors in the exchange energy. 
These errors are absent in QMC calculations. They can also be viewed as arising from the nonanalytic behavior of 
the HF structure factor at k = 0. We have constructed an accurate correction for these errors in HF theory. 

For 2D systems the leading-order finite-size errors (using both the Ewald and MPC interactions) are caused by the 
slow convergence of the XC hole and the neglect of long-ranged correlations in the KE. The errors in the energy per 
particle scale as 0(iV -5 / 4 ), suggesting that this form should be assumed in the extrapolation to infinite system size. 

If the single-particle finite-size error is to be removed by extrapolation rather than twist averaging then the HF 
exchange energy should not be included in the extrapolation; just the KE. Furthermore, an estimate of the effective 
mass should be included in the extrapolation. 

Tests at realistic system sizes show that time-step bias in DMC results does not get significantly worse as the system 
size is increased. Population control bias does get worse, but only slowly. 
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APPENDIX A: FINITE-SIZE ERRORS IN HF THEORY 

For the 3D HEG, the HF exchange hole is±£ 

2 



(Al) 

> 1/3 

where in the last line we have retained only the dominant nonoscillatory term at large separation, fcp CT = (67r 2 iV -/f2j 
is the Fermi wave vector for particles of spin a, and N a is the number of particles of spin a. The hole has a slowly 
decaying tail that falls off as 1/r 4 , so there is a missing contribution to the exchange energy in a finite simulation cell. 
The interaction of each electron with its exchange hole should be 1/r (as enforced inside the simulation cell when the 
MPC interaction is used). So the missing contribution to the HF interaction energy is approximately 

a<> » 

1 JRa r 

= " i^r^' 3 . <A2) 

where Rq is the radius of a sphere of volume il. This gives a finite-size error in the HF energy per particle that falls 
off slowly as iV" 2 / 3 . This error will also be present in the Ewald energy. In addition to this missing contribution, 
there are errors arising from the fact that the part of the exchange hole that would lie outside the simulation cell if 
the system were infinite is distorted by being compressed back into the simulation cell to satisfy the sum rule. The 
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charge of the missing tail is approximately 

Q = [°° ^r 2 p x {r) dr = — (^) V N 2 J 3 — = # ( —Y V ^ /3 - (A3) 

If we assume that this missing charge is uniformly distributed inside a sphere of radius we must subtract its 
unwanted contribution to the exchange energy, giving another correction 

( 2 ) _ _WQ (J_Y /3 ^ N 2/3 (M) 

AVhf - W Q ~ 2wr s UttJV J ^ N °- (A4j 

(7 

(Other approximations, such as assuming Q to increase linearly within Rq may be more accurate.) The total correction 
to the exchange energy (either Ewald or MPC) obtained within this real-space procedure is 

AV HF = AV^ + AV$ = -L (JL) V3 £ NT. (A5) 

The result of applying this correction to the HF Ewald exchange energy is shown in Fig. [6l along with the result of 
applying the correction of Eq. (|4ip . Both work well, although the correction of Eq. (f4"T|) is better. 



APPENDIX B: EQUIVALENCE OF THE MPC AND XC CORRECTION 

Consider a system of cubic symmetry. The difference between the MPC and Ewald XC energies is: 



(Vmpc) - (Vew) = y J /0 xc (r) f ~ - [vs(r) ~ vm}^ dr 



- t L ft = (r) (I- 3 + •■■)*■ <B1) 

were we have used the expansion of the Ewald interaction from Eq. (|19p . Assuming that p xc (r) is well localized within 
the simulation cell, we can replace /? X c(r) by S\ oc (y) — 5(r) and extend the range of integration to infinity to obtain 



(Vmpc)-(Vew) = S / r 2 S loc (r)dr 

>J" Jr<oo 



= -f^V£S loc (k)| k=0 . (B2) 
Since 5(k) = r]k 2 + 0(k 4 ) in a cubic system, this reproduces Eq. ([25)) : 

(Vmpc) - (VW) = + ■ • ■ ■ (B3) 

The use of the first-order AVe w correction may therefore be regarded as a first-order approximation to the MPC, in 
which the leading term in the small-r expansion of the difference between 1/r and v_e(r) — vm is taken into account 
but higher order terms are neglected. 
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